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We study numerically the disorder-induced localization-delocalization phase transitions that oc- 
cur for mass and spring constant disorder in a three-dimensional cubic lattice with harmonic cou- 
plings. We show that, while the phase diagrams exhibit regions of stable and unstable waves, the 
universality of the transitions is the same for mass and spring constant disorder throughout all 
the phase boundaries. The combined value for the critical exponent of the localization lengths of 
v — 1.550^o'oi7 confirms the agreement with the universality class of the standard electronic An- 
derson model of localization. We further support our investigation with studies of the density of 
states, the participation numbers and wave function statistics. 



PACS numbers: 63.50.-x, 63.20.D-, 63.20.Pw 

I. INTRODUCTION 

The disorder-induced metal-insulator transition (MIT) 
and the concept of Anderson localization 1 have been 
studied extensively for over 50 years. Most of the atten- 
tion was focused on electronic systems and their trans- 
port propertied^' — indeed the acronym MIT itself sug- 
gests this. However, localization physics is of course 
much broader than just electrons in solid state devices 
and encompasses the whole realm of waves — quantum 
and classical — and their interference due to random 
scattering events. Recently, the interest in localization 
has been rekindled by its beautiful realization in cold 
atom systems.— 7 Similarly, localization of classical waves 
has received new impetus from spatially resolved studies 
in elastic, vibrational systems. 8 

Theoretical work on the localization properties of har- 
monic solids has received somewhat less attention over 
the years. In our opinion, this could be due to (i) a 
general expectation that the vibrational problem only 
mimics the electronic one and (ii) the one clear feature 
when this is not the case — the so-called "boson peak" 
(BP) 9 10 — up to this date remains to be understood 
fully. In a recent paper p] we have shown that expecta- 
tion (i) is only partially true: the phase diagrams, even 
for just a simple cubic harmonic lattice of masses and 
springs, exhibit several intriguing features for both the 
purely mass and the purely spring constant disordered 
cases. A similarly distinguishing characteristic of vibra- 
tional localization is the fact that the zero frequency, i.e. 
uo = mode that corresponds to global translational in- 
variance . ca nnot be localized regardless of the amount of 
disorder.^ The aforementioned BP corresponds to the 
appearance of a low-frequency enhancement of the den- 
sity of states g(uo) with respect to Debye's g(uS) oc uo 2 
law.™ Mogt 

previous investigations of the localization 
properties of disordered vibrational modes agre e that the 
modes near and above the BP are extended) 9 * 13 * 14 * i.e. 
cjbp <C co> c , where ujbp denotes the BP frequency (peak 



of g(uo)/uu 2 ) and uu c the boundary between extended and 
localized states. It has been argued before via eigenvalue 
statistics that the states with ujbp < uo < uj c are governed 
by random-matrix statistics of the Gaussian orthogonal 
ensemble (GOE).™ 



In this paper, we present a detailed study of the vibra- 
tional localization and transport properties throughout 
the previously obtained phase diagrams of a cubic har- 
monic lattice system with either random mass or random 
spring constant disorder. Using large matrix diagonal- 
ization techniques, we investigate the behaviour of the 
vibrational density of states (VDOS) as well as the par- 
ticipation numbers and wave function statistics of the 
vibrational eigenstates. This complements earlier stud- 
ies of participation ratios JSE3 level- spacing statistical^ 
and multifractal properties.^ In particular we demon- 
strate that the disorder-effected states below uj c exhibit 
a modified Porter-Thomas statistics of the wave func- 
tions, which is close to the one from the GOE ensem- 
ble. In addition, we present results from a high-precision 
transfer matrix method (TMM) and a finite-size scaling 
(FSS) analysis which allow us to corroborate the phase 
diagrams and calculate the universality class of the mo- 
bility edges across all of the phase diagram. Our results 
have relevance in the related problem of instantaneous 
normal modes in glasses and supercooled liquids^^ as 
well as acoustic metamaterialsP^^ Here we just note 
that in both these classes of materials, there exist exci- 
tations which can be related to the existence of states in 
what is formally part of the temporally decaying, neg- 
ative uj 2 region of the phase diagrams shown in Fig. [I] 
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FIG. 1. (Color online) Phase diagrams for (a) mass (Am vs. uj 2 ) and (b) spring constant (A A: vs. uj 2 ) disorder. Grey and 

(dark) brown shadings are the regions of extended and localized states, respectively. The white regions are inaccessible and their 

perimeter denote the band edges. The horizontal dash-dotted line indicates the border between stable and unstable regions, 

the dotted lines denote uj 2 — and 6. Thick black crosses, white circles, open blue circles and closed green circles indicate (i) 

the three high-precision transition values obtained from FSS, (ii) the other lower- accuracy transition points, (iii) the reduced 

VDOS peak locations from numerics and (iv) the reduced VDOS peak locations computed via CPA, respectively, in both (a) 

and (b). The yellow labels @-© show the positions of the states in Figs. a) -(c) and Figs. |3ja)-(c). The dashed blue lines 

indicate the range of deviations Eq. |7| from universal wave function statistics, cp. also the insets of Figs. [6] and [7] Additionally 

in (a) the (very dark) marroon shading denotes the critical region obtained from the transformation of the electronic phase 
diagram.™ 



II. SCALAR MODEL OF LATTICE 
VIBRATIONS 

A. The clean case 

We shall consider masses arranged on a simple cubic 
lattice and connected by harmonic forces. With Uj de- 
noting the deviation from the lattice equilibrium position 
fj = (x, y, z)j of a certain mass rnj at given x, y and z 
lattice coordinates, we can write the classical equations 
of motion as 

/ k x \ / u x \ 
m jtij = - k v \ \ u v ' (*) 

all neighbours n 

where k xi k yi k z and u x , u yi u z denote the spring con- 
stants and displacements in y and z direction for each 
nearest neighbour n, respectively. Often the components 
of the spring constant are categorised into central and 
n on- central terms, central when acting along the dimen- 
sion of their subscript, e.g., k x along the x-direction and 
non- central otherwise. We can reduce the computational 
complexity of the problem by assuming that central and 
non-central force constants are identical. This turns all 
force constant matrices into scalars. After this reduction 
the three dimensions of the system are decoupled into 
three identical independent problems and solving any one 
solves the ful l syst em. This "scalar" model, or "isotropic 
Born model" j * 32 * can be written in its stationary form 



as 

- uj 2 mjUj = y^kji(ui - Uj), (2) 
i 

where uj is the frequency of vibration and Uj(t) = Uje lUJt . 
In matrix notation, we have an eigensystem with eigen- 
values —uj 2 , 

- uj 2 U = M _1 KU, (3) 

where M _1 K is called the dynamical matrix and, due to 
infinitesimal translational symmetry^ always obeys the 
sum rule ^^(M _1 K) J 7 = 0. In the clean case, we have 
that all masses are equal to a constant m and all spring 
constants are k. With these definitions, the frequencies 
range from to the largest possible frequency o;^ ax = 
12k /m and uj 2 will always be given in units of [k/m]. 



B. The disordered case 

We are interested in introducing disorder into the sys- 
tem. From ([3]), it is clear that this can be done (i) 
by allowing the masses to vary such that rnj G [fn — 
Am/2, m + Am/2] and (ii) by having random spring con- 
stants kji G [k — Ak/2,k + Afc/2]. For simplicity, we 
will use the uniform mass and spring constant distribu- 
tions with m = k = 1 and restrict our investigation to 
the two cases of either pure mass or pure spring con- 
stant disorder. Note that this choice sets the units as 



3 



well. The classical problem presented in Eq. ([!]), partic- 
ularly its stationary form (J2|, is very similar to the tight- 
binding Schrodinger equation for the three-dimensional 
Anderson model of localization 1 at energy E such that 
(E—€j)ipj = — J2 t tjiipi, where the / summation is over all 
nearest neighbours and ej and tji denote the onsite and 
hopping energies, respectively!^ For the mass-disordered 
model with fluctuating masses rrij one can obtain the 
transformation relations 



E 6- 



6j(E) ^ uj 2 mj = (6 - E)rrij . (4) 



As shown in Ref. 11, we can then reuse many of the 
results for the Anderson model and infer the phase di- 
agrams of localization-delocalization transitions for the 
vibrational mass-disorder model. In Fig. [TJa), we show 
the estimated mobility edges for the case of pure vi- 
brational mass disorder based on transforming the re- 
lated estimates of the mobility edges in the Anderson 
modelP^HU The phase diagrams for the vibrational case 
are intriguing in many respects.^ First of all (i) there is 
clear evidence for delocalization-localization transitions 
due to disorder. Next, (ii) the strong disorder limits of 
1 2 Am | > m, with the possibility of negative masses, or 
|2Afc| > fc, with similarly possible negative spring con- 
stants, give rise to locally unstable regions (although 
globally stable) corresponding to negative uo 2 solutions. 
Such modes are known in liquids as unstable instanta- 
neous normal modes and are related to the relaxation 
dynamics of the liquids (iii) The separation of ex- 
tended and localized states continues into these regions 
and hence do the transitions and (iv) there is a re-entrant 
behaviour for uo > and Am (Afc) < 2. These extraordi- 
nary mobility edges and hence the phase diagra ms hav e 
been confirmed by direct high-precision numerics! 11 * 37 * 



III. LOCALIZATION PROPERTIES OF 
EIGENSTATES 

A. Numerical diagonalization 

Let us start our investigation of Q by looking at some 
typical eigenstates obtained by exact diagonalization. In 
particular, we are using a combination of the iterative nu- 
merical eigensystem packages ARPACK^and PardisoP^ 
We find this combination to be most effective when deal- 
ing with both the unsymmetric and the symmetric cases 
of pure mass and spring disorder, respectively.^ 

In Fig. [2j we show eigenstates for the pure mass dis- 
order case corresponding to three eigenfrequencies which 
lie in regions which according to the phase diagram (Fig. 
[]Ja)) should be extended, close to the mobility edge and 
localized. We see from Fig.|2]that these characterisations 
reflect the apparent nature of these vibrational states. 
For Fig. [2|a), the local amplitude of vibrations at each 
site is roughly of similar magnitude throughout the sys- 
tem, whereas for Fig. [2|c), the vibrations are confined 
to a small region in the cube. Figure |5|b) displays the 



characteristic properties of a critical wave function at the 
Anderson mobility edge.™ 

For the pure spring disorder case as in Fig. [3j we see 
that the vibrations for the three shown frequency values 
may also be classified into extended, critical and local- 
ized classes. This classification indeed agrees with the 
computed phase diagram as shown in Fig. [TJb) for the 
pure spring disorder case. However, we also see that 
the character of the states seem subtly different from 
the pure mass disorder ones. The vibrations seem to 
be more around certain vibration centres and radiate 
outward roughly symmetrically from these centres!^ Al- 
though not the topic of the present investigation, we em- 
phasise that this should make the multifractal analysis of 
such states very informative, in particular its comparison 
with the recently proposed symmetry of the multifractal 
spectrum! 17 1 41 1 43 1 



B. Vibrational density of states 

In order to numerically obtain the VDOS, the compu- 
tation of all states is required. The iterative methods 
applied in section III A| are then no longer efficient and 
we employ a standard LaPackP^ dense matrix routine 
(DGEEV). 

We have calculated the VDOS g(uj 2 ) = g(uj)/2uj for 
disorders Am, Afc = 0.5, 1, 1.5, 2, 2.5, 3, 4, 5, 6, 7, 8, 9 and 
10 for cubes with volume L 3 = 15 3 for 50 disorder con- 
figurations. This results in roughly 170, 000 Uj's for each 
disorder. In Figs. [4ja) and pj^b) we show the results 
for g(uj 2 ) as a functions of w for all mass and spring 
constant disorder magnitudes respectively. We find for 
both types of disorder that the van Hove singularities 
in the VDOS become smeared out upon increasing the 
disorder. In addition, there are the usual low-frequency 
peaks corresponding to standing waves in the simula- 
tion box. These pea ks indicate the presence of plane- 
wave-like states J=^"=3 We perform analytical calculations 
of the VDOS using the coherent-potential approximation 
(CPA, see Appendix A) .^Except for the standing- wave 
peaks (which are absent in the L — »• oo CPA calcula- 
tions) there is very good agreement between the analyt- 
ical and numerical results as can be seen from Figs. Qa) 
andgb). Using the CPA one can easily evaluate the 



maxima of the "reduced VDOS" g(uj)/uj 2 = 2g(uj 2 )/uj 
("boson peaks"). For small disorder (Am < 1, Afc < 1) 
these peaks are identical with the transverse van Hove 
singularities, located at uo 2 =4. For larger disorder the 
BPs become disorder-dominated and no longer reflect 
the underlying lattice symmetry. This can be (and has 
been) checked by CPA calculations using a Debye Green's 
function Gq(z) = gD(ty/(z — A) with A = a; 2 , 
gD{u 2 ) = 3cj/2cj|). In these calculations the BP posi- 
tions for Am > 1, Afc > 1 coincide with those of the 
lattice calculations. It has been shown in Ref. [9] that the 
BP separates a nearly plain wave regime from a regime 
where disorder is dominant (random-matrix regime). We 
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FIG. 2. (Color online) Schematic representation of amplitude distributions \uj \ obtained from exact diagonalization for system 
of length L 3 = 70 3 for mass disorder Am = 4 and frequencies (a) uo 2 = 3, (b) uJ 2 = 4.5 and (c) u 2 = 6. All sites with 
u(fj) I L 3 ^2j u (?j) > 1 are shown as small cubes and those with black edges have u(fj) / L 3 ^2j u(fj) > a/1000. The color scale 
distinguishes between different slices of the system along the axis into the page. 




FIG. 3. (Color online) Schematic representation of amplitude distributions \uj \ as in Fig. [2] obtained for spring disorder Ak = 1 
and frequencies (a) uo 2 — 12, (b) uo 2 — 12.5 and (c) uo 2 — 13.03. The size of the cubes and their colour is chosen as in Fig. ^1 
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FIG. 4. (Color online) Vibrational density of states g(uo 2 ) as 
a function of frequency cj 2 and various disorders of (a) Am 
and (b) Ak. The blue and red lines in the base denote the 
trajectories of the localization-delocalization transition and 
the band edges, respectively, see Fig.[l] The thin dashed lines 
for Am, Ak = 5, 10 are results from CPA calculations (see 
Appendix [A] for details) . The thick dashed line is the clean 
simple cubic density of states and identical in both cases. 



find from analysing our VDOS data that this is also the 
case for our model systems. 

However the scenario for mass and spring-constant dis- 
order is very different. In the spring-constant disorder 
case the BP, and with it, the range of nearly plane waves 
goes continuously towards zero near Ak = 2.5. In CPA 
there are no states with uj 2 < below this value. In the 
mass disorder case the BPs and correspondingly the low- 
frequency range of nearly plane waves extend towards 
Am — » oo. This can be easily understood by the transfor- 
mation rule Q, which states that the mass fluctuations 
are suppressed by a factor uj 2 . Therefore for uj 2 — » 
there are always plane waves in the infinite- volume sys- 
tem, which are converted to standing waves at finite vol- 
ume. 

It is interesting to note that in the mass disorder case 
a peak on the negative uj 2 side develops for high values 
of Am near the uj 2 < mobility edge. On the positive 
uo 2 side both the peak in g(co> 2 ), the BP and the mobil- 
ity edge approach each other with increasing Am. This 
confirms that there is no proportionality between co>bp 
and uj c as postulated in Ref. [lOj The absence of such a 



simple relationship was also already discussed in Refs. 
and|50l 



C. Participation numbers 

The participation number Pl(cj u ) is a measure of the 
number of sites in the lattice that are contributing to the 
vibrational excitation of the nth vibrational eigenstate 
u\ (n) , U2 (n) , . . . , u L 3 (n) . It can be denned aJ^ 



PZ\<»n) 



3 



uAn) 



(5) 



in analogy with the electronic case. We emphasise that 
the normalisation Y^j u 2 (n) = 1, automatically observed 
for electronic eigenstates by the Born rule, has to be en- 
forced for the vibrational case for consistency in the com- 
parison between different eigenstates.^ A fully extended 
vibration will lead to Pz,(a; n ) — 1 whereas a vibration 
localized at a single site corresponds to Pl{^u) = l/L 3 
and hence in the limit L — » oo. 

We average the participation numbers in discrete fre- 
quency intervals over 50 disorder realizations and plot 
them for each disorder at L = 15 in Figs. [5ja) and |5|b) 
for mass and spring disorder, respectively. We find that 
the transition from delocalized to localized behaviour as 
found in section |IV A| does not lead to a clear crossing of 
Pl for system sizes L = 5, 10 and 15. We expect to see 
such a crossing only when going to much larger system 
sizes and upon increasing the number of disorder sam- 
ples. Thus our results show the difficulties associated 
with the use of participation numbers in studying the 
present transition in agreement with a recent attempt by 
Monthus et alP^ 1 

In general, the results for Pl nevertheless confirm the 
phase diagrams presented in Figs. [IJa) and [ljb) as the 
extended regions of the phase diagrams are matched with 
states of higher participation. The VDOS results of Fig. 
[4] are also confirmed qualitatively as extended states usu- 
ally lead to higher Pl values than localized ones. In par- 
ticular, we note the emergence of finite Pl values in the 
negative u 2 regime for large mass disorder as well as the 
pronounced tail in the same frequency regime for strong 
spring constant disorder. 



D. Vibrational Eigenstate Statistics 

Disordered quantum systems exhibit irregular fluctu- 
ations of eigenfunctions, which ca n be studied from the 
statistics of the local amplitudes! 52 * 53 * In the universal 
regime (of mostly weak disorder), random matrix the- 
ory can classify these fluctuations into universality classes 
such as the Porter-Thomas distribution^] Q f the GOE.^ 
Upon increasing the disorder, corrections to GOE have 
been studied which we expect t o see present also in the 
case of our vibrational disorder! 56 * 57 * We determine the 
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/ and fg OE as 



FIG. 5. (Color online) Participation numbers Pi 5 as a func- 
tion of uj 2 for various (a) mass disorders Am and (b) spring 
constant disorders Ak, averaged over 50 disorder realisations. 
The dotted grey and red lines in the base denote the phase 
boundaries and the band edges respectively as in Fig. [I] 



distribution function 



L 3 



J25(v-\ Uj (n)\*L 3 )5(^ 



(6) 



where A is the mean level-spacing, ( ) denotes an average 
over disorder realisations and the vibrational eigenvectors 
are normalised so that (\u j(n)\ 2 ) = L~ 3 . For all disor- 
ders mentioned in section IIIB we calculate f(v) from 
over two million amplitudes at L 3 = 70 3 for frequencies 
throughout the phase diagram at intervals of Suj 2 = 0.5 
and plot them for mass and spring disorder in Figs. 
[6|a)-[7ja), respectively. We include the Wigner estimate 
from random matrix theory, /qq E = exp(— v/2) / yj2wvW^ 
For exponentially localized states, one finds f^(v) ~ 
In 2 (c 2 L 3 1 v) j 'v with c a dis orde r dependent constant and 
£ is the localisation length OSED We also include the the- 
oretical result for a maximally localized scenario, where 
£ = 1 in Figs. [6ja)-[7^a). We see in Fig. J^a) that the 
curves for increasing mass disorder increasingly depart 
from GOE, whereas for spring constant disorder in Fig. 
[TJa) there is an abrupt departure from GOE when the 
localization-delocalization transition is crossed. In Figs. 
[6^b) and[7jb) we plot the relative difference Sf between 



f(v) 



Sf(v) = 



(7) 



GOE 



and include the analytical estimate of departure from 
GOE as derived for the electronic Anderson model 56 



Sf(v) 



3v v 2 

T + T 



(8) 



where A is a constant related to the diffusion in the sys- 
tem. We see that for small frequencies the analytical 
estimate is very well suited to our data and we show that 
for the mass disorder case a value of A\ = 0.0545 has 
a good fit for Sf of uj 2 = 2 and similarly A\ = 0.0315 
has a good fit for Sf of uj 2 = 2.5 in the spring constant 
disorder case. For higher frequencies this fit continues 
in the spring constant disorder case, where for a value 
A2 = 0.0545 we have a good agreement with Sf of uj 2 = 6. 
This is not the case in the mass disorder case where the 
minimum values of Sf shift from v = 3 and as an illus- 
tration we show that for A2 = 0.195 the difference Sf 
fits the uj 2 = 3.5 results only for small v but very quickly 
deviates for increasing v. 

Upon further increasing cj 2 , we see that there is again 
a region where the agreement with /q^ e becomes better. 
This behaviour has not previously been observed (neither 
in the electronic case nor in calculations on vibrational 
modes). 

In the inserts of Figs.[6^b) and[7jb) we have plotted the 
minima of Sf(v) as a function of frequency for different 
values of the disorder parameters Am and Ak. As stated 
above, these functions exhibit a minimum corresponding 
to a maximum deviation of the eigenstate fluctuations 
from the GOE behaviour. We have marked the positions 
of these minima in the phase diagrams in Fig. [T] and find 
that they coincide with the values of the BP frequencies. 
Obviously both the disorder-modified plane waves (uj < 
cjbp) as well as the random-matrix states (uj > cjbp) 
obey the GOE statistics rather well, whereas the states 
at the cross-over (i.e. the states with uj = cjbp) have a 
maximum deviation from GOE. Of course, approaching 
the mobility edge the GOE behaviour disappears. 



IV. LOCALIZATION PROPERTIES OF 
TRANSPORT STATES 

A. TMM results and the phase diagrams for mass 
and spring disorder 

We have performed TMM calculations at Am, Ak = 
0.2, 0.4, . . . , 2 (see Appendix [B] for details). In addition 
to these disorders, more are required to verify the phase 
boundary obtained for the pure mass disorder case from 
direct transformation of the electronic potential disor- 

A small selection 



dered phase boundary in section II B 



of additional disorders is chosen as Am = 2.2,4,6,9. In 
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FIG. 6. (Color online) For disorder Am = 1 (a) f(v) for a 
range of frequencies as labelled in the figure, with the pure 
GOE result indicated by the dashed line, and the highly lo- 
calized result indicated by the dot-dashed line, (b) Sf(v) for a 
range of frequencies where every 5 th data point has a symbol. 
GOE here corresponds to the Sf = line whereas the vertical 
v — 3 line is the minimum position of the analytical Sf. The 
dot-dashed and dotted lines indicate Sf of ^ with constants 
A\ — 0.0545 and A2 = 0.195, respectively. Inset: the value of 
the minimum of 5f plotted as a function of uj 2 for different 
mass disorders. The green dots correspond to the positions 
of the BP obtained via CPA calculations and the dashed blue 
lines estimate the position of the minima and their width is 
equal to the separation of the minimum function data points, 
all shown in the phase diagram in Fig. [TJa). 



the pure spring disordered case a larger additional list is 
required as a phase boundary is yet to be established. An 
adequate resolution is achieved with additional disorders 
of Ak = 2.5, 3, 4, 4.5, 5, 6, 7, 8, 9, 10. The average of the 
mass and spring constant disorder (m and k) has been 
kept fixed at 1 for all cases. For every disorder value, the 
reduced localization length, Am has been calculated for 
a range of frequencies and system widths M = 6, 8, 10 
and 12 to an accuracy of 0.1% of the variance. 

In Figs. [8(a]j|8(b)| we show the resulting disorder and 




dependencies for 2 of the 6 representative mass/spring 
disorder regions. At all disorder magnitudes for both 
spring constant and mass disorder, these figures reveal 



FIG. 7. (Color online) For disorder A A; = 1 (a) f(v) for a 
range of frequencies as labelled in the figure with the pure 
GOE and fully localized behaviours as indicated in Fig.[6|a). 
Panel (b) shows Sf{v) for a range of frequencies where every 
5 th data point has a symbol. GOE here corresponds to the 
Sf — line whereas the vertical v = 3 line is the minimum 
position of the analytical Sf. The dot-dashed and dotted 
lines indicate Sf of |8| with constants A\ — 0.0315 and A2 = 
0.0545, respectively. Inset: as in Fig. [6] also plotted in the 
corresponding phase diagram Fig. [TJb). 



clear transitions from extended behaviour, with increas- 
ing Am values for increasing M, to localized behaviour, 
where Am decreases when M increases. We also see 
in these figures frequency regions where Am remains 
roughly constant upon changing M. Such regions are 
in the vicinity of a change from derealization to local- 
ization and hence Figs. |8(a)-[8(bJ indicate the existence 
of a delocalization-localization transition. We roughly 
estimate the transition regions by the frequency value at 
which the values of Am for the largest and the second 
largest system size cross (M = 10, 12). Then we obtain 
a similarly rough estimate of the error of this estimate 
from the difference with respect to the frequency value 
which we obtain when we take the crossing point between 
the largest and smallest system sizes (M = 6, 12). These 
estimates are the basis of the phase diagrams in Fig. [I] 

In the spring constant disorder case we ne ed t o pay 
special attention to the k~+ x term in equation ( |B2| ) as at 
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(a) rir = 2, n ri =3, n % = 1, m r = 2, m 2 = 



(b) n ro = 2, n ri = 2, 



FIG. 8. (Color online) Reduced localization lengths Am plotted as function of uj 2 for various system sizes as indicated by 
different symbols. Panel (a) shows mass disorder for Am = 1.2 while (b) is for spring constant disorder of Ak = 7. The lines 
in each plot show the fits obtained from FSS, the orders of the expansion are given below each figure (see Appendix [C]) . The 
vertical dotted line represents the estimated values of uj 2 with orange shading indicating the error obtained from Monte Carlo 
analysis (in Tab. l|. Error bars are only shown for the largest and smallest system size, as in all cases they are within symbol 
size. The insets display the obtained scaling function when the irrelevant components have been subtracted. 



disorders Ak > 2, the disorder distribution contains val- 
ues close to zero which when applied in the k~ +1 can 
dramatically increase a single site amplitude dwarfing 
surrounding amplitudes. We apply a cut-off whereby if 
\k x +i\ < 10 _4 /c the value is rejected and another ran- 
domly chosen. We re-estimate all transition frequencies 
of previously mentioned disorders and find that the new 
estimates are identical within the previous error bars and 
therefore keep the estimates obtained with unaltered dis- 
tributions. 

We plot these estimates of the critical frequencies in 
the phase diagrams of Figs. [TJa) and[TJb). As we can 
see, for the pure mass disorder case, Fig. [IJa) very well 
reproduces the estimated phase diagram obtained from 
comparison with the electronic phase diagram in the An- 
derson model.^ Most interestingly, the small pocket of 
extended states in the complex frequency spectrum of the 
mass disorder phase diagram is clearly identified by the 
two transitions from localized to delocalized and back to 
localized at Am = 9. 

For the pure spring constant disorder, we see that in 
the region < uj 2 < 12, all states remain extended up to 
the largest considered spring constant disorder Afc = 10. 
This is similar to the electronic case with pure hopping 
disordei^M^ where even very strong hopping disorder 
does not lead to complete localization close to E = 

We find that both for m ass and spring constant disor- 
der, the uj 2 = mode^ULL^ES remains extended regardless 
of the disorder strength. This is in agreement with pre- 
vious studies in one- and two-dimensional systems.^ We 
also observe for both mass and spring constant disor- 
der very strong shifts of the crossing points of Am when 
changing M. This is to be expected since we are effec- 
tively dealing with transition regions in the vicinity of the 



tails of the VDOS (cp. Fig. [4| and hence the systematic 
size changes are also strongly influenced by non-universal 
changes in the VDOS. This is similar to the situation 
for the electronic case where the transition at the mo- 
bility edges for E ^ is known to be more difficult to 
studyPfe] 



B. FSS estimates for the critical parameters 

In order to obtain more reliable estimates for the tran- 
sition point uj 2 as well as to ascertain the existence of a di- 



vergent correlation length £(cj) ex \uj 



at uj n with 



critical exponent ^, we need to proceed to the M — >• oo 
limit. This we do, as in the electronic case, via an FSS 
procedure (see Appendix |c| for details) We perform 
the FSS analysis on the raw data of reduced localiza- 
tion lengths Am as functions of uj 2 as well as uj (with 
£(cj) oc \uo — uo c \~ v ). While the latter seems more natural 
in the context of vibrations, we emphasise that the for- 
mer is more convenient when comparing to the electronic 
case where uj 2 is related to the energy. 11 

For both pure mass and pure spring disorder, we con- 
centrate on 3 disorder values each, choosing those from 
the 3 different domains of the phase diagrams of Figs. 
[TJa) andQb), namely (i) uj 2 > 0; Am, Ak < 2, (ii) 
uj 2 > 0; Am, A/c > 2 and (iii) uj 2 < 0. For these 6 
points, we compute additional high-precision data for 
M = 14, 16, 18 and 20. The additional Am values fo r two 
o f thes e 6 transitions have also been shown in Figs. 8(aJ 



- |8(b)[ We then apply the FSS procedure of appendix [C 
and hence obtain precise estimates of the critical param- 
eters and transition frequencies uj 2 of a vibrating solid 
in the thermodynamic limit. These uj 2 values have also 
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been indicated in the phase diagrams as in Figs.[]Ja) and 
[TJb)). In Tab. [i] we show the results for the high-precision 
FSS analysis. We find that in all cases, a consistent, ro- 
bust and stable fit with quality-of-fit parameter T q larger 
than 0.1 can be identified. In particular, the FSS for lj 
as well as uo 2 gives consistent results. 



A weighted average of the critical exponent for the es- 
timates in Tab [i] is v — 1.550^o'oi7- This is in excellent 
agreement with previous numerical studies of the Ander- 
son model for electron localization which have found the 
critical exponent v = 1.57 ± 0.02pH^In the vibrational 
model, no previous high-precision results are available. 
With an accuracy of 2% in the raw TMM data for spring 
disorder Ak = 1.8, Akita and OhtsukJ^p rev k)usly found 
a critical exponent of v ~ 1.2 ±0.2. Recently, Monthus 
and Garel 51 assumed v = 1.57 and showed that their 
participation ratio data for high disorder collapsed onto a 
scaling function. All these results for model Q are there- 
fore consistent with the orthogonal universality class of 
the Anderson model P 



V. CONCLUSIONS 



Making contact with possible experimental systems, 
we note that the transitions are at rather high frequen- 
cies. The Debye temperatures = ftwjj/kB of, e.g., 
Si and Ge — ca ndidate materials for milli-Kelvin cool- 
ing device^E^ whose study got us interested in this 
research — are Qd = 645K and 374K, respectively. 
Assuming that the upper band edge of the clean case 
can be approximated by the respective Debye frequen- 
cies uj d = 1.34 x 10 13 Hz and 7.79 x 10 12 Hz, respectively, 
we see from the phase diagrams that the transition fre- 
quencies remain quite high. Localization of vibrations 
for these systems in the stable regime appears only pos- 
sible for frequencies in or above the far infrared frequency 
spectrum, particularly for spring constant disorder. The 
transition for very large mass disorder does tend towards 
smaller u 2 values, but these mass disorders are already 
deep in the unstable regime Am > 2. This is of course 
dramatically different from the electronic situation where 
a disorder of 16.55 is known to localize all states in a 
simple cubic system with band width 12 (in units of hop- 
ping strength) P We note that the unstable regions of the 
phase diagrams for Am, Ak > 2 with possibly negative 
masses and spring constants are now recognised to be of 
considerable interest fo r aco ustic and disordered meta- 
material applications mmMM}! 

ere our identification of 
regions of extended states should prove useful. 



In the preceding sections, we have established the ex- 
istence and universality of the localization-delocalization 
transitions for vibrational excitations in a simple har- 
monic solid at various values of frequency and mass or 
spring constant disorder. While the model itself is simple, 
the resulting phase diagrams are not and exhibit intrigu- 
ing features. In particular, there are regions of localized 
and extended unstable modes with transitions between 
them that belong to the same universality class as in the 
stable regimes. Namely, the universality class of the 3D 
electronic Anderson metal-insulator transition!^ Our re- 
sults show that the FSS scaling works both when using 
the uj scaling, most natural from a vibrational point of 
view, as well as the uo 2 scaling, motivated by the elec- 
tronic analogue. The peak in the VDOS as shown in Fig. 
H] seems identifiable as a continuation of the van Hove 
singularity at low — mass or spring constant — disor- 
der. The peak is not visible in the participation ratio 
data, but its signature can be seen again in the wave 
function statistics. Whether it can truly be called a bo- 
son peak, although it does of course appears as such in 
g(uj)/uj 2 plots, remains undetermined at present.^ The 



wave function statistics of section HID and the plots of 
critical vibrational amplitudes in Fig. [2] and [3] also reveal 
subtle differences between mass and spring disorder. A 
more in-depth analysis of the multifractal properties and 
scaling properties of the generalised participation ratio 
at the transition might be very useful. However, we note 
that previous studies in fluids^ and elastic beadd^ have 
found good agreement with t he mu ltifractal spectrum ob- 
tained for the electronic case. 41 69 
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Appendix A: Coherent potential approximation 



As an estimate of the VDOS calculations of section 



IIIB| we compute the VDOS using the coherent potential 
approximation (CPA).^In the spring constant disorder 
case we introduce a frequency dependant force constant 
("self energy") T(z) and determine its contribution self- 
consistently using the scattering matrix formalism, 9 



T(z) - h 



1 _ [r^-kij^i-zGjz)] 
1 3r(*) 



= 0, 



(Al) 



where z = uj 2 + ^0 + is the regularised complex frequency. 
The local Green function of the effective medium is 



(A2) 



and Go is the Green function for the clean system!^ The 
averaged VDOS is then given as 



(g(aj 2 )) = --3m[G(z)}. 

7T 



(A3) 
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Am M 



1.2 8 - 20 [12.15, 13.1] 2 3 12 

4.0 8 - 20 [3.75, 4.25] 3 2 110 

9.0 8 - 20 [-1.65, -1.5] 2 3 12 

1.2 8 - 20 [3.485, 3.62] 2 3 12 

4.0 8 - 20 [1.936, 2.062] 3 2 110 

9.0 8 - 20 [-1.284, -1.225] 2 3 110 



12.681 



+0.056 
0.034 
a -,o 4 +0.024 
4: -- LO4: -0.020 
1 £ 9 q+0.018 
±.UZO_q q 37 



+0.14 
0.09 



o r R1 +0.008 
O.OD1_ 005 

5+0.006 
-0.005 
1 970+0. 006 
-l.Z <Q_Q 014 



2.033Z 



1.57 

r 7 +0.06 

1 "° ' -0.08 
+0.41 
0.18 
+0.15 
0.09 

-| rr+0.07 

1 - oo -0.08 
-+0.44 
-0.17 



1.56 
1.57 



1.561 



165112 
572^ 

164t 3 3 * 
573t 6 6 * 
155_ 33 



165 0.84 
574 0.99 

154 0.87 
165 0.84 
573 0.99 

155 0.83 



Ak M 



1.0 


10 - 


20 




[12.48, 12.6] 


3 


1 


1 


1 


1 




19 ^97+0-003 
' -0.004 


1 co + 0.05 
l.OO_ .04 


132^ 


132 0.62 


10.0 


6 - 


16 




[18.8, 20.3] 


1 


3 


1 


2 







IQ 74Q+ - 043 

±y. i ^_o.038 


1 +0.08 
1 ' O1 -0.08 




176 0.84 


7.0 


8 - 


20 




[-3.5, -2.75] 


2 


2 


1 


1 







q qOr+0.070 

°- ozo -0.115 


1 5Q+0-23 


162111 


162 0.51 


1.0 


10 - 


20 


[3.529, 3.55] 




3 


3 


1 


1 


2 


3.5401°:°°! 




1 47+0.15 

-0.05 


15711J 


156 0.49 


10.0 


6 - 


16 


[4.335, 4.506] 




2 


3 


1 


2 





4.4411°;™ 




c-9+0.15 
1 ' OZ -0.53 


199l£ 


199 0.87 


7.0 


8 - 


20 


[-1.87, -1.66] 




2 


2 


1 


1 





±.O^U_ 0- o33 




1 60 +0 - 21 
i.ou_ 19 


162lg 


162 0.79 



TABLE I. Values of critical parameters cj c , &c an d v for pure mass (top) and pure spring constant (bottom) disorder computed 



from FSS performed in the given M and cj, uj ranges and with the orders of the expansion (CI ) given by n ro , n ri , ni, m r and 
mi. The minimised % 2 value, the degrees of freedom \i and the resulting goodness-of-fit parameter T q are also shown for each 
fit. The errors correspond to non-symmetric 95% confidence intervals (see Appendix |C|) . 



In the mass disordered case we use the transformation 
rule Q to map the problem to an Anderson problem with 
fluctuating local energies e$ and then use the conventional 
single-site CPAP^The self energy E(z) with z = E + i0+ 
is given by setting the following CPA scattering matrix 
equal to zero: 



E(s) 



l-(e io -E(*))G [z-£(*) 







(A4) 



The single-site CPA problem is known to exhibit rather 
unstable iteration properties. We obtained a good itera- 
tion performance using the following iteration method^, 
which is equivalent to the CPA condition (A4). 



E^ n+1 \z) =S(")( Z ) + 



(t) 



(«) 



(n) 



1+ (t)( n )G [z-T,( n )(z)] 

;„-£<»>(*) 



l-[e i0 -S(»)( 2 )]G [^-S(»)(z)] 



, (A5) 



(A6) 



where n is the iteration count. The average density of 
states is then calculated from the Green's function as 



(9(E)) 



1 



3m {G [*-£(*)]}. 



(A7) 



The results for both disorders are shown in Fig. [4] as 
thin dashed lines next to the numerical VDOS. We find 
good agreement between CPA results and the numerical 
calculations for both weak and strong disorder and in the 
stable (uj 2 > 0) and unstable (uj 2 < 0) spectral regions. 



Appendix B: The transfer- matrix approach 

The transfer-matrix method (TMM) allows for a very 
memory efficient way to iteratively calculate the decay 



length Am of vibrations in a quasi-one dimensional bar 
with cross section M x M for lengths L ^> M. Equation 
Q has to be rearranged into a form where the amplitude 
of vibration of a site in layer x + 1 — when x is chosen 
as the direction of transfer — is calculated solely from 
parameters of sites in previous layers x and x — 1, 



1 



^x-\-l,y,z 



kx-\-l,y 

kx — l,y,z 
kx+l,y,z 



[(uj 2 m X} y }Z + K\\)u Xi y, z - h x 



Ux-1,- 



(Bl) 



Here h 

k. 



x,y,z-\-l^x,y,z-\-l ~r ^x^y^z — l^x^y^z — l \ 

'■x,y+i,zV>x,y+i,z + ^x,y-\,z^x,y-\,z denotes the collection 
of in-plane contributions to the final amplitude, k a \\ = 

kx,y,z-\-l H~ kx,y,z — 1 ~^~kx,y— l,z H~ k x + l,y,z l,y,z 

and we have changed back to the explicit notation 
such that uj = u Xj y jZ for fj = (x,y,z)j. Similarly, 
kji k x+ljVjZ for ri = {x + l,y,z)i. With U x 
(^,1,1,^,1,2,^,2,1, • • • ,u XiM ,m), we can define U x , U x+1 
and U x -\ as vectors containing the amplitudes of the 
constituent sites in l ayers x, x + 1 and x — 1, respec- 
tively. Equation (Bl) can now be expressed in standard 



transfer-matrix form 



U x +i 




u x 





[(c^m^+fcan)!-^] 
k x +i 



k x +i 





±1 



u x 

Ux-l 



(B2) 

where is a M x M matrix containing all in-layer con- 
tributions, and 1 are the zero and unit matrices, re- 
spectively. 

Formally, the transfer matrix is used to 'trans- 
fer' vibrational amplitudes U from one slice to the 
next and repeated multiplication of this gives the global 
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transfer matrix tl = Yl^i^x- The limiting matrix 

i 

T = lim^^oo (tl,t\)j 2L exists^ and has eigenvalues e ±7i , 

i = 1,...,M. The inverse of these Lyapunov expo- 
nents are estimates of decay /localization lengths and 
the physically relevant largest vibrational decay length 
is Am(^ 2 ) = 1/min^ ^(cj 2 )]. The reduced (dimension- 
less) decay length may then be calculated as Am(w 2 ) = 
A M (w 2 )/M. 



Appendix C: Finite-size scaling 

The FSS includes two types of corrections to scaling, 
namely, those which account for the nonlinearities of the 
Am, Ak dependence of the scaling variables (relevant 
scaling) and for the mentioned shift of the point at which 
the Am(w 2 ) curves cross (irrelevant scaling). The start- 
ing point for the FSS in terms of uo 2 is the scaling ansatz 

A M (Lo 2 ) = f(xrMi, X iMy), (CI) 

where Xr and Xi are the relevant and irrelevant scaling 
variables, respectively. The function Am is then Tay- 
lor expanded up to the order rii and we have Am = 
12n=o XiM ny f n (xrMi^ from where we obtain a series 
of functions f n which are in turn Taylor expanded up to 
an order n r such that f n (xrM^J = Ylk=o a nkXr M " • 



Nonlinearities are taken into account by expanding both 
Xi and Xr in terms of the dimensionless frequency w = 
(u 2 -uj 2 )/uj 2 such that Xr(w) = YlZ=i h m™ m , XiO) = 
Ylm=o c m w7n where the orders of the expansions are m r 
and rrii. For a more rigorous analysis we hard-code the 
zero-th and first order of the irrelevant expansion and 
Taylor expand each appearance of f n separately!^ 

The expansions of the fit functions and the fit are per- 
formed numerically up to the orders n ro , n ri , , rrii and 
m r . Each individual data set can be best suited to a 
particular expansion, the general rule being that the or- 
ders of expansion should be kept as low as possible while 
giving the best fit to the data, and minimising the esti- 
mated standard errors for the critical parameters u 2 and 
v. We check for stability of the fit by individually in- 
creasing each expansion parameter by one and checking 
to see that the obtained parameters remain within the 
95% confidence intervals of the original fit. 

The confidence intervals are then recomputed through 
a Monte Carlo analysis.^ We obtain a perfect data se- 
ries from the fit with the previously computed expansion. 
We next vary each data point according to a Gaussian 
distribution with the right standard deviation. With this 
synthetic data, we then repeat the FSS fit to obtain new 
estimates of the critical parameter. We repeat this oper- 
ation 5000 times and compute the distribution function 
for each critical parameter. We then estimate the true 
errors from these histograms by taking as errors those 
values at which 2.5% of the distribution are below or 
above bulk. 
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